Particle in a Box¶
Particle-in-a-box: $$ \hat{H} = \frac{-\hbar^2}{2m}\frac{\partial^2}{\partial x^2} + V(x) $$
where $0<x<L = V$, and $x=\infty$ everywhere else. Solve to find:
$$ \begin{align} \psi(x) &= \sqrt{\frac{2}{L}}\sin(\frac{n\pi}{L} x) \\ E_n = \frac{n^2\hbar^2\pi^2}{2mL^2} \end{align} $$
If we include the time-dependent phase factor:
$$ \Psi(x,t) = \psi(x)e^{-iEt/\hbar} $$
Constants:
- reduced Planck's constant: $(\hbar) = 1.05457182*10^{-34} \; \text{m}^2 \text{kg s}^{-1}$
- mass of electron: $m_e = 9.1093837*10^{-31} \; \text{kg}$
- angstrom: $\text{\AA} = 10^{-10} \; \text{m}$
In [ ]:
import matplotlib
# Enable interactive plot
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import numpy as np
A = 1e-10 # m; one angstrom
me = 9.1093837e-31 # kg; mass of electron
hbar = 1.05457182e-34 # m^2 kg / s
n_grid = 100
# L = 10*A
#
L = 1
hbar = 1
x = np.linspace(0, L, n_grid)
delta_x = x[1]-x[0]
def energy(n, m=me, L=L):
return (n*hbar*np.pi)**2 / (2*m*L**2)
def norm(x, p):
pass
def psi(x, n, t=0, m=me, L=L):
E = energy(n, m, L)
return np.sqrt(2/L)*np.sin((np.pi*n*x)/L)*np.exp((-1j*E*t)/hbar)
Plot¶
Now let's plot the wavefunction
In [ ]:
import plotly.express as px
import pandas as pd
# t = np.linspace(0, 1e-15, 100)
t = np.linspace(0, 10, 100)
data = []
# m = me
m = 1
ki = 2 # this is the eigenstate
print(" Energy is: ", energy(1, m=m, L=L))
for ti in t:
nt = 0
for xi in x:
p = psi(xi, ki, t=ti, m=m)
data.append({'x':xi, 'Real':p.real, 'Imag':p.imag, 'Dens':(p.conj()*p).real, 't':ti})
df = pd.DataFrame(data)
ymax = np.max(np.abs(df.Dens))
fig = px.scatter(df, x='x', y=['Imag', 'Real', 'Dens'], animation_frame='t', range_y=[-ymax,ymax])
print(np.sum(df[df.t==0]["Dens"])*delta_x)
fig.show()
Energy is: 4.934802200544679 1.0
Non-stationary state¶
Now let's plot the time evolution of the system in a superposition of different energy eigenstates:
In [ ]:
data = []
m = 1
k = [1,2,3,4,5,6,7,8,9,10]
k = [1,2]
print(" Energy is: ", energy(1, m=m, L=L))
for ti in t:
nt = 0
for xi in x:
p = 0
for ki in k:
p += psi(xi, ki, t=ti, m=m)
p = p / np.sqrt(len(k))
data.append({'x':xi, 'Real':p.real, 'Imag':p.imag, 'Dens':(p.conj()*p).real, 't':ti})
df = pd.DataFrame(data)
ymax = np.max(np.abs(df.Dens))
fig = px.scatter(df, x='x', y=['Imag', 'Real', 'Dens'], animation_frame='t', range_y=[-ymax,ymax])
fig.show()
Energy is: 4.934802200544679
In [ ]:
fig = px.line_3d(df, x='x', y='Imag', z='Real', labels={'x':'x', 'y':'Real', 'z':'Imag'}, animation_frame="t", range_y=[-ymax,ymax], range_z=[-ymax,ymax])
fig.show()